The goal was to determine if particles on the electron microscope images follow spiral or concetric arrangement.
import matplotlib as mpl
from matplotlib import pyplot as plt
import numpy as np
import cv2
import os
import copy
from scipy.optimize import curve_fit
import nanoparticles as nano
import importlib
importlib.reload(nano)
<module 'nanoparticles' from 'h:\\My Drive\\Projekty\\Small projects\\spiral\\nanoparticles.py'>
file = 'W14_J1_3_47 [1].tif'
img = cv2.imread('data/' + file)
mask = cv2.imread('data/' + file.replace('.', '_mask.'))
mask = cv2.cvtColor(mask, cv2.COLOR_RGB2GRAY)
centers = cv2.imread('data/' + file.replace('.', '_centers.'))
img = cv2.flip(img, 1)
mask = cv2.flip(mask, 1)
centers = cv2.flip(centers, 1)
imgGray = cv2.cvtColor(img.copy(), cv2.COLOR_RGB2GRAY)
scale, img_scale = nano.find_scalebar(
imgGray.copy(), 200, showscale=True) # scale is in nm/px
print(f'Scale is {scale:.4} nm/px')
nano.plot_images([img_scale, mask, centers], names=['Raw Image (+scalebar)',
'Manually selected regions', 'Manually selected centers of spirals'], size=14)
Scale is 0.3035 nm/px
def get_spirals_centers(centers):
centersGray = cv2.cvtColor(centers, cv2.COLOR_RGB2GRAY)
ret, thrCen = cv2.threshold(255-centersGray, 5, 255, cv2.THRESH_BINARY)
contours, hierarchy = cv2.findContours(
thrCen, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
spirals = []
# colored = cv2.cvtColor(imgGray,cv2.COLOR_GRAY2RGB)
for c in contours:
# compute the center of the contour
M = cv2.moments(c)
cX = int(M["m10"] / M["m00"])
cY = int(M["m01"] / M["m00"])
spirals.append([cX, cY])
return spirals
spirals = get_spirals_centers(centers)
fewSpirals = [nano.get_aoi(img, spiral) for spiral in spirals[:5]]
nano.plot_images(fewSpirals, size=8, rowsplit=5)
blur = cv2.medianBlur(imgGray, 9)
# use adaptive gaussian thresholding
blocksize = 251
thr = cv2.adaptiveThreshold(
blur, 255, cv2.ADAPTIVE_THRESH_GAUSSIAN_C, cv2.THRESH_BINARY, blocksize, 2)
thr = np.where(thr == 255, 0, 255).astype(np.uint8)
# erode particles
kernel = np.ones((5, 5), np.uint8)
thrEroded = cv2.erode(thr.copy(), kernel, iterations=1)
# select particles
dMin, dMax = 1, 7 # min and max diameter of particles in nm
gCon, bCon = nano.find_contours(
thrEroded, scale, areaMin=np.pi*(dMin/2)**2, areaMax=np.pi*(dMax/2)**2, circMin=0.5)
rgbColors = cv2.cvtColor(imgGray,cv2.COLOR_GRAY2RGB)
rgbColors = cv2.drawContours(rgbColors, gCon, -1, (0,255,0), 1)
rgbColors = cv2.drawContours(rgbColors, bCon, -1, (255,0,0), 1)
particles = nano.get_centers(gCon)
names = ['Blurred', 'Threshold', 'Eroded', 'Particles']
nano.plot_images([nano.get_aoi(image, spirals[1], size = 500) for image in [blur, thr, thrEroded, rgbColors]], names = names, size = 12)
# filter particles that are outside mask
spiralPrtcs = []
for p in particles:
if mask[p[1], p[0]] != 0:
spiralPrtcs.append(p)
spiralPrtcs = np.array(spiralPrtcs)
# assign to spiral based on distance to centers
particlesAssigned = [[] for _ in range(len(spirals))]
for p in spiralPrtcs:
distances = []
for s in spirals:
distances.append(nano.point_distance(p, s))
minD = np.argmin(np.array(distances))
particlesAssigned[minD].append(p)
plt.imshow(imgGray, cmap = 'gray')
for prt in particlesAssigned:
plt.scatter(np.array(prt)[:, 0], np.array(prt)[:, 1], s = 0.3)
plt.scatter(np.array(spirals)[:, 0], np.array(spirals)[:, 1], c = 'black', s = 15)
for i, spiral in enumerate(spirals):
x, y = spiral[0], spiral[1]
plt.text(x, y, str(i), color = 'black')
rgbColors = cv2.cvtColor(imgGray,cv2.COLOR_GRAY2RGB)
for n in range(len(spirals)):
color = nano.random_color()
for particle in particlesAssigned[n]:
rgbColors = cv2.circle(rgbColors, particle, 5, color, -1)
nano.plot_images([nano.get_aoi(rgbColors, spiral) for spiral in spirals[:4]], size = 12)
# add polar coordinates
particlesDetails = []
for i in range(len(spirals)):
spiral = spirals[i]
particles = []
for x, y in particlesAssigned[i]:
dist, alfa = nano.cart2pol(x-spiral[0], y-spiral[1])
particles.append([x, y, dist, alfa])
particlesDetails.append(np.array(particles))
We are looking for rows of particles that follow along the curve around the center of spiral. We first select one particle (that is closes to center) and then find next particle like this:
How alfa was calculated
maxDist = 3.5/scale
idealSpiralDist = 9/scale
drs = []
# grouping particles into rows that follow curve around the center of spiral
def group_particles(particlesDetails):
particlesGrouped = []
for i in range(len(particlesDetails)):
group = []
particles = copy.copy(particlesDetails[i])
particles = np.array(particles)
particles = list(particles[particles[:, 2].argsort()]) # sort by distance to center
prt = particles[0] # get first particle
row = [prt] # creates initial row
del particles[0]
direction = 0 # direction (clockwise or anti) should be unifrom within row of particles
while len(particles) > 0:
x, y, r, theta = row[-1]
# alfa is the expect angle of ideal neighbouring particle at that particular distance to center (r)
alfa = 2*np.arcsin(idealSpiralDist/(2*r))
# if direction is selected
if direction != 0:
x2, y2 = nano.pol2cart(r, theta+alfa*direction)
dsts = [nano.point_distance((p[0] - spirals[i][0], p[1] - spirals[i][1]), (x2, y2)) for p in particles]
minD = np.argmin(dsts)
if dsts[minD] < maxDist:
row.append(particles[minD])
del particles[minD]
else:
if len(row) > 1: group.append(row) # only take rows with more than 1 particle
row = [particles[0]] # starts new row
del particles[0]
direction = 0
else: # find direction (check on which side best particle match is closest to ideal)
# positive direction
x2P, y2P = nano.pol2cart(r, theta+alfa)
dstsP = [nano.point_distance((p[0] - spirals[i][0], p[1] - spirals[i][1]), (x2P, y2P)) for p in particles]
minDP = np.argmin(dstsP)
# negative direction
x2N, y2N = nano.pol2cart(r, theta-alfa)
dstsN = [nano.point_distance((p[0] - spirals[i][0], p[1] - spirals[i][1]), (x2N, y2N)) for p in particles]
minDN = np.argmin(dstsN)
# select
if min(dstsN[minDN], dstsP[minDP]) < maxDist:
if dstsN[minDN] < dstsP[minDP]:
row.append(particles[minDN])
del particles[minDN]
direction = -1
else:
row.append(particles[minDP])
del particles[minDP]
direction = 1
drs.append(direction)
else:
row = [particles[0]]
del particles[0]
direction = 0
particlesGrouped.append(group)
return particlesGrouped
particlesGrouped = group_particles(particlesDetails)
C:\Users\szust\AppData\Local\Temp\ipykernel_23676\666208268.py:21: RuntimeWarning: invalid value encountered in arcsin alfa = 2*np.arcsin(idealSpiralDist/(2*r))
# just drawing
work = copy.copy(img)
for n, spiral in enumerate(particlesGrouped):
work = cv2.circle(work, spirals[n], 10, (255, 255, 255), -1)
for group in spiral:
color = nano.random_color()
for i in range(0, len(group)-1):
p1 = np.array(group[i][:2])
p2 = np.array(group[i+1][:2])
work = cv2.line(work, p1.astype(int), p2.astype(int), color, 2)
work = cv2.circle(work, p1.astype(int), 4, color, -1)
work = cv2.circle(work, p2.astype(int), 4, color, -1)
nano.plot_images([nano.get_aoi(work, spiral) for spiral in spirals[:4]], size = 12)
def linear_func(x, a, b):
return a*x + b
def fix_angles_inrow(angles):
'''
add/substract 2*np.pi from list of angles after going through 0/2*np.pi (idicated when theta difference is > np.pi)
'''
newAngles = [angles[0]]
adder = 0
for i in range(1, len(angles)):
if angles[i-1] - angles[i] > np.pi:
adder += 2*np.pi
if angles[i] - angles[i-1] > np.pi:
adder -= 2*np.pi
newAngles.append(angles[i] + adder)
return np.array(newAngles)
def calc_and_plot(rows, work, center):
fig, axs = plt.subplots(ncols = 2, nrows = 1, figsize = (12, 6*1))
factorRng = 5 # to determine color
coefs, thetasSelected, distancesSelected = [], [], []
for row in rows:
row = np.array(row)
if len(row) > 2 and min(row[:, 2])*scale > 10: # take only rows that are not too close to center
distances = row[:, 2]
thetas = row[:, 3]
thetas += np.pi # to not have negative theta (does not matter really)
thetas = fix_angles_inrow(thetas)
# fitting to linear equation to get "a" coeficient (y = ax + b)
popt, pcov = curve_fit(linear_func, thetas, distances, bounds=(-np.inf, np.inf))
a = popt[0]*scale
coefs.append(a)
meanCoef = np.mean(coefs)
# draw this row on chart
c = mpl.cm.get_cmap('coolwarm')((factorRng+a)/(2*factorRng))
axs[0].text(np.mean(thetas), np.mean(distances)*scale, f'{a*scale:.1f}', fontsize = 7)
axs[0].scatter(thetas, distances*scale, s = 6, color = c)
axs[0].plot(thetas, distances*scale, color = c, lw = 1.7)
axs[0].set_xlabel('Theta [rad]')
axs[0].set_ylabel('Distance [nm]')
axs[0].set_title('Fitted a (dist = a * theta + b)')
# draw this row on image
ccv = (np.array(c[:3])*255).astype(int) # color used by cv2
ccv = tuple(map(int, ccv))
for ip in range(0, len(row)-1):
p1 = np.array(row[ip][:2])
p2 = np.array(row[ip+1][:2])
work = cv2.line(work, p1.astype(int), p2.astype(int), ccv, 5)
work = cv2.circle(work, p1.astype(int), int(2/scale), ccv, -2)
work = cv2.circle(work, p2.astype(int), int(2/scale), ccv, -2)
# draw colored dot at center of spiral
c = mpl.cm.get_cmap('coolwarm')((factorRng+meanCoef)/(2*factorRng))
ccv = (np.array(c[:3])*255).astype(int) # color used by cv2
ccv = tuple(map(int, ccv))
work = cv2.circle(work, center, int(4/scale), (0,0,0), -2)
work = cv2.circle(work, center, int(3/scale), (ccv), -2)
axs[1].imshow(nano.get_aoi(work, center))
axs[1].set_title('Corresponding image')
thetasSelected.append(thetas)
distancesSelected.append(distances)
plt.close()
return fig, [coefs, thetasSelected, distancesSelected]
fig, result = calc_and_plot(particlesGrouped[2], copy.copy(img), spirals[2])
fig
fig, result = calc_and_plot(particlesGrouped[6], copy.copy(img), spirals[6])
fig
fig, result = calc_and_plot(particlesGrouped[8], copy.copy(img), spirals[8])
fig
work = copy.copy(img)
results = []
for i in range(len(spirals)):
fig, result = calc_and_plot(particlesGrouped[i], work, spirals[i])
results.append(result)
nano.plot_images([work], size = 30)
fig, axs = plt.subplots(ncols = 2, nrows = 1, figsize = (12, 6*1))
factorRng = 5
allCoefs = []
for result in results:
coefs, thetasSelected, distancesSelected = copy.copy(result)
for coef, thetas, distances in zip(coefs, thetasSelected, distancesSelected):
c = mpl.cm.get_cmap('coolwarm')((factorRng+coef)/(2*factorRng))
firstPoint = np.argmin(thetas)
thetas -= thetas[firstPoint]
distances -= distances[firstPoint]
axs[0].scatter(thetas, distances*scale, s = 6, color = c)
axs[0].plot(thetas, distances*scale, color = c, lw = 1.7)
axs[0].set_xlabel('Theta [rad]')
axs[0].set_ylabel('ΔR [nm]')
axs[0].set_title('Fitted coefficient (R = a * theta + b)')
allCoefs.append(coef)
_ = axs[1].hist(allCoefs, bins = 40)
axs[1].set_title('Coefficient histogram')
axs[1].set_xlabel('Coefficient')
axs[1].set_ylabel('Count')
Text(0, 0.5, 'Count')
CONCLUSION: Depending on the specific spiral, particles can follow right (a > 0, red color) or left (a < 0, blue color) spiral arrangement or be close to concetric (a = 0).